///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/lbfgs.h
///
/// \brief   header file for the LBFGS method
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef FORMIC_LBFGS_HEADER
#define FORMIC_LBFGS_HEADER

#include<vector>
#include<numeric>
#include<cassert>
#include<algorithm>
#include<cmath>

#include<boost/math/special_functions/fpclassify.hpp>

#include<formic/mpi/interface.h>
#include<formic/output/output.h>
#include<formic/exception.h>
#include<formic/lapack/interface.h>
#include<formic/numeric/numeric.h>

namespace formic {

//  template <class Manager>
//  void lbfgs(const int n,
//             double & f,
//             double * const x,
//             double * const g,
//             const int max_iter,
//             const int max_hist,
//             Manager & man) {
//
//    std::vector<double> rho(max_hist, 0.0);
//    std::vector<double> alpha(max_hist, 0.0);
//    std::vector<double> beta(max_hist, 0.0);
//
//    std::vector<double> old_x_vec(n, 0.0);
//    std::vector<double> old_g_vec(n, 0.0);
//    double * const old_x = &old_x_vec.at(0);
//    double * const old_g = &old_g_vec.at(0);
//
//    for (int iter = 0; iter < max_iter; iter++) {
//
//      // save old parameter and gradient vectors
//      if ( iter > 0 ) {
//        formic::xcopy(n, x, 1, old_x, 1);
//        formic::xcopy(n, g, 1, old_g, 1);
//      }
//
//      // get the value
//      man.get_val(n, x, f);
//
//      // get the gradient
//      man.get_grd(n, x, g);
//
//      // compute differences
//      if ( iter > 0 ) {
//        formic::xscal(n, -1.0, old_x, 1);
//        formic::xscal(n, -1.0, old_g, 1);
//      }
//
//      // compute the rho vector
//      for (int i = 1; i < pmtr_hist.size() && i <= userinp.get<int>("bfgs_hist_len"); i++) {
//        //rho.at(i) += formic::xdotc(pmtr_hist.at(i-1)->size(), &pmtr_hist.at(i-1)->at(0), 1, &grad_hist.at(i-1)->at(0), 1);
//        //rho.at(i) -= formic::xdotc(pmtr_hist.at(i-1)->size(), &pmtr_hist.at(i-0)->at(0), 1, &grad_hist.at(i-1)->at(0), 1);
//        //rho.at(i) -= formic::xdotc(pmtr_hist.at(i-1)->size(), &pmtr_hist.at(i-1)->at(0), 1, &grad_hist.at(i-0)->at(0), 1);
//        //rho.at(i) += formic::xdotc(pmtr_hist.at(i-1)->size(), &pmtr_hist.at(i-0)->at(0), 1, &grad_hist.at(i-0)->at(0), 1);
//        //rho.at(i) = 1.0 / rho.at(i);
//        formic::xcopy(pdiff.size(),       &pmtr_hist.at(i-1)->at(0), 1, &pdiff.at(0), 1);
//        formic::xaxpy(pdiff.size(), -1.0, &pmtr_hist.at(i-0)->at(0), 1, &pdiff.at(0), 1);
//        formic::xcopy(gdiff.size(),       &grad_hist.at(i-1)->at(0), 1, &gdiff.at(0), 1);
//        formic::xaxpy(gdiff.size(), -1.0, &grad_hist.at(i-0)->at(0), 1, &gdiff.at(0), 1);
//        rho.at(i) = 1.0 / formic::xdotc(pdiff.size(), &pdiff.at(0), 1, &gdiff.at(0), 1);
//        if (formic::mpi::rank() == 0)
//          formic::of << boost::format("rho[%2i] = %.2e") % i % rho.at(i) << std::endl;
//      }
//
//      // prepare space for the alpha and beta vectors
//
//      // compute the update
//      std::vector<double> update(**grad_hist.begin());
//      for (int i = 1; i < pmtr_hist.size() && i <= userinp.get<int>("bfgs_hist_len"); i++) {
//        //alpha.at(i) += rho.at(i) * formic::xdotc(update.size(), &pmtr_hist.at(i-1)->at(0), 1, &update.at(0), 1);
//        //alpha.at(i) -= rho.at(i) * formic::xdotc(update.size(), &pmtr_hist.at(i-0)->at(0), 1, &update.at(0), 1);
//        //formic::xaxpy(update.size(), -alpha.at(i), &grad_hist.at(i-1)->at(0), 1, &update.at(0), 1);
//        //formic::xaxpy(update.size(),  alpha.at(i), &grad_hist.at(i-0)->at(0), 1, &update.at(0), 1);
//        formic::xcopy(pdiff.size(),       &pmtr_hist.at(i-1)->at(0), 1, &pdiff.at(0), 1);
//        formic::xaxpy(pdiff.size(), -1.0, &pmtr_hist.at(i-0)->at(0), 1, &pdiff.at(0), 1);
//        formic::xcopy(gdiff.size(),       &grad_hist.at(i-1)->at(0), 1, &gdiff.at(0), 1);
//        formic::xaxpy(gdiff.size(), -1.0, &grad_hist.at(i-0)->at(0), 1, &gdiff.at(0), 1);
//        alpha.at(i) = rho.at(i) * formic::xdotc(update.size(), &pdiff.at(0), 1, &update.at(0), 1);
//        formic::xaxpy(update.size(), -alpha.at(i), &gdiff.at(0), 1, &update.at(0), 1);
//      }
//      for (int i = std::min(int(pmtr_hist.size() - 1), userinp.get<int>("bfgs_hist_len")); i > 0; i--) {
//        //beta.at(i) += rho.at(i) * formic::xdotc(update.size(), &grad_hist.at(i-1)->at(0), 1, &update.at(0), 1);
//        //beta.at(i) -= rho.at(i) * formic::xdotc(update.size(), &grad_hist.at(i-0)->at(0), 1, &update.at(0), 1);
//        //formic::xaxpy(update.size(), alpha.at(i) - beta.at(i), &pmtr_hist.at(i-1)->at(0), 1, &update.at(0), 1);
//        //formic::xaxpy(update.size(), beta.at(i) - alpha.at(i), &pmtr_hist.at(i-0)->at(0), 1, &update.at(0), 1);
//        formic::xcopy(pdiff.size(),       &pmtr_hist.at(i-1)->at(0), 1, &pdiff.at(0), 1);
//        formic::xaxpy(pdiff.size(), -1.0, &pmtr_hist.at(i-0)->at(0), 1, &pdiff.at(0), 1);
//        formic::xcopy(gdiff.size(),       &grad_hist.at(i-1)->at(0), 1, &gdiff.at(0), 1);
//        formic::xaxpy(gdiff.size(), -1.0, &grad_hist.at(i-0)->at(0), 1, &gdiff.at(0), 1);
//        beta.at(i) = rho.at(i) * formic::xdotc(update.size(), &gdiff.at(0), 1, &update.at(0), 1);
//        formic::xaxpy(update.size(), alpha.at(i)-beta.at(i), &pdiff.at(0), 1, &update.at(0), 1);
//      }
//      formic::xscal(update.size(), -1.0, &update.at(0), 1);
//
//      // convert the update to complex form if necessary
//      std::vector<S> cplx_update(co->n_der(), formic::zero(S()));
//      local_funcs::real_to_cplx(co->n_der(), &update.at(0), &cplx_update.at(0));
//
//      // perform a line search along the update direction and remember the optimum step size
//      formic::mpi::bcast(_step_size);
//      double step_arg = _step_size;
//      formic::fqmc::line_search(userinp, burn_len, samp_len, start_config, start_rng, &cplx_update.at(0), userinp.get<double>("update_max_sdev"), old_co, co, step_arg);
//      if (step_arg > 0.0)
//        _step_size = step_arg;
//
//    }
//
//  }

  template<class W>
  inline void quadratic_line_search(W & worker,
                                    const int n,
                                    const double f0,
                                    const double * const x0,
                                    const double * const p,
                                    double * const x1,
                                    double & step_arg,
                                    const bool bcast_step_size = true) {

    //// start a timer
    //formic::start_timer("lbgfs loose line search");

    // get the initial value
    const double init_energy = f0;

    // print initial value
    if (formic::mpi::rank() == 0) {
      formic::of << boost::format("for a step size of %16.8e") % 0.0;
      formic::of << boost::format("      val = %20.12f") % init_energy << std::endl << std::endl;
    }

    // initialize line search variables
    std::vector<double> step_sizes(1, 0.0);
    std::vector<double> step_energies(1, init_energy);
    double best_step_size = 0.0;
    double best_step_energy = init_energy;
    int best_step = 0;
    int step_i = 0;
    int n_refine = 0;

    for (int snum = 1; snum <= 20; snum++) {

      double new_step_size = -1.0;
      bool to_break = false;

      // choose what to do next
      if ( step_sizes.size() < 2 ) {
        new_step_size = step_arg;
        if (formic::mpi::rank() == 0)
          formic::of << boost::format("trying first step...") << std::endl << std::endl;
      } else if ( step_sizes.size() < 3 ) {
        new_step_size = step_sizes.at(step_i) * ( best_step_energy == init_energy ? 0.45 : 1.3 );
        if (formic::mpi::rank() == 0)
          formic::of << boost::format("trying second step...") << std::endl << std::endl;
      } else if ( best_step == 0 ) {
        new_step_size = step_sizes.at(1) * 0.5;
        if (formic::mpi::rank() == 0)
          formic::of << boost::format("so far all steps are too long, so we will try a shorter step") << std::endl << std::endl;
      } else if ( best_step == step_sizes.size()-1
                  && *step_energies.rbegin() < *(step_energies.rbegin()+1)
                     + ( *step_sizes.rbegin() - *(step_sizes.rbegin()+1) ) * ( *(step_energies.rbegin()+1) - *(step_energies.rbegin()+2) ) / ( *(step_sizes.rbegin()+1) - *(step_sizes.rbegin()+2) )
                ) {
        // ez < ey + ( z - y ) * ( ey - ex ) / ( y - x )
        if (formic::mpi::rank() == 0)
          formic::of << boost::format("stopping line search due to negative curvature") << std::endl << std::endl;
        to_break = true;
      } else {
        if ( best_step == step_sizes.size()-1 )
          new_step_size = formic::porabola_min_max(&step_sizes.at(best_step-2), &step_energies.at(best_step-2));
        else
          new_step_size = formic::porabola_min_max(&step_sizes.at(best_step-1), &step_energies.at(best_step-1));
        new_step_size = std::min(new_step_size, 2.0 * step_sizes.at(best_step));
        if (formic::mpi::rank() == 0)
          formic::of << boost::format("choose final step based on a quadratic fit") << std::endl << std::endl;
      }

      formic::mpi::bcast(to_break);
      if (to_break)
        break;

      // if requested, broadcast the step size to all processes
      if (bcast_step_size)
        formic::mpi::bcast(new_step_size);

      // prepare the next step
      for (step_i = 0; step_i < step_sizes.size(); step_i++) {
        if ( new_step_size < step_sizes.at(step_i) )
          break;
      }
      step_sizes.insert(step_sizes.begin()+step_i, new_step_size);
      step_energies.insert(step_energies.begin()+step_i, 1.0e100);
      if ( step_i <= best_step )
        best_step++;

      // compute new value
      if (formic::mpi::rank() == 0) {
        formic::xcopy(n, x0, 1, x1, 1);
        formic::xaxpy(n, step_sizes.at(step_i), p, 1, x1, 1);
      }
      worker.get_value(n, x1, step_energies.at(step_i));
      if (formic::mpi::rank() == 0) {
        formic::of << boost::format("for a step size of %16.8e") % step_sizes.at(step_i);
        formic::of << boost::format("      val = %20.12f") % step_energies.at(step_i) << std::endl << std::endl;
      }

      // if this was the best step size so far, save it
      if ( step_energies.at(step_i) < best_step_energy ) {
        best_step = step_i;
        best_step_size = step_sizes.at(step_i);
        best_step_energy = step_energies.at(step_i);
      }

      formic::mpi::bcast(step_sizes);
      formic::mpi::bcast(step_energies);
      formic::mpi::bcast(best_step);
      formic::mpi::bcast(best_step_size);
      formic::mpi::bcast(best_step_energy);

      // check if we are finished
      if ( step_sizes.size() > 3 && best_step != 0 )
        break;

    }

    if (formic::mpi::rank() == 0) {

      // if a zero step was best, do not apply any update
      if ( best_step == 0 ) {

        formic::xcopy(n, x0, 1, x1, 1);

        formic::of << "The line search failed to find a step that lowers the energy.  No update was applied." << std::endl << std::endl;

      // otherwise, apply the update corresponding to the best step size
      } else {

        formic::xcopy(n, x0, 1, x1, 1);
        formic::xaxpy(n, best_step_size, p, 1, x1, 1);

        if (formic::mpi::rank() == 0)
          formic::of << boost::format("The line search's best step size was %12.6f with a value of %20.12f") % best_step_size % best_step_energy
                     << std::endl
                     << boost::format("The corresponding update has been applied.")
                     << std::endl
                     << std::endl;

      }

    }

    //// stop the timer
    //formic::stop_timer("lbgfs loose line search");

    // return the step size used
    step_arg = best_step_size;

  }

  template<class W>
  inline void lbfgs_loose_line_search(W & worker,
                                      const int n,
                                      const double * const x0,
                                      const double * const p,
                                      double * const x1,
                                      double & step_arg,
                                      const bool bcast_step_size = true) {

    //// start a timer
    //formic::start_timer("lbgfs loose line search");

    // get the initial value
    double init_energy;
    worker.get_value(n, x0, init_energy);

    // print initial value
    if (formic::mpi::rank() == 0) {
      formic::of << boost::format("for a step size of %16.8e") % 0.0;
      formic::of << boost::format("      val = %20.12f") % init_energy << std::endl << std::endl;
    }

    // initialize line search variables
    std::vector<double> step_sizes(1, 0.0);
    std::vector<double> step_energies(1, init_energy);
    double best_step_size = 0.0;
    double best_step_energy = init_energy;
    int best_step = 0;
    int step_i = 0;
    int n_refine = 0;

    for (int snum = 0; snum < 30; snum++) {

      double new_step_size = -1.0;

      // choose what to do next
      if ( step_sizes.size() < 2 ) {
        new_step_size = step_arg;
        if (formic::mpi::rank() == 0)
          formic::of << boost::format("trying first step...") << std::endl << std::endl;
      } else if ( step_sizes.size() < 3 ) {
        new_step_size = step_sizes.at(step_i) * ( best_step_energy == init_energy ? 0.5 : 1.3 );
        if (formic::mpi::rank() == 0)
          formic::of << boost::format("trying second step...") << std::endl << std::endl;
      } else if ( best_step == 0 ) {
        new_step_size = step_sizes.at(1) * 0.5;
        if (formic::mpi::rank() == 0)
          formic::of << boost::format("so far all steps are too long, so we will try a shorter step") << std::endl << std::endl;
      } else if ( best_step == step_sizes.size()-1 ) {
        new_step_size = best_step_size * 1.3;
        if (formic::mpi::rank() == 0)
          formic::of << boost::format("best step so far was longest step, so we will try an even longer step") << std::endl << std::endl;
      } else {
        n_refine++;
        new_step_size = formic::porabola_min_max(&step_sizes.at(best_step-1), &step_energies.at(best_step-1));
        // make sure the new step size is well behaved
        if ( new_step_size < 0.0 || boost::math::isnan(new_step_size) ) {
          break; // stop if we have found a negative number or a nan
        } else if ( boost::math::isinf(new_step_size) || new_step_size > 1.4 * (*step_sizes.rbegin()) ) {
          new_step_size = 1.4 * (*step_sizes.rbegin());
        } else if ( new_step_size < 0.6 * step_sizes.at(1) ) {
          new_step_size = 0.6 * step_sizes.at(1);
        }
        if (formic::mpi::rank() == 0)
          formic::of << boost::format("refining the local minima found during the line search") << std::endl << std::endl;
      }

      // if requested, broadcast the step size to all processes
      if (bcast_step_size)
        formic::mpi::bcast(new_step_size);

      // prepare the next step
      for (step_i = 0; step_i < step_sizes.size(); step_i++) {
        if ( new_step_size < step_sizes.at(step_i) )
          break;
      }
      step_sizes.insert(step_sizes.begin()+step_i, new_step_size);
      step_energies.insert(step_energies.begin()+step_i, 1.0e100);
      if ( step_i <= best_step )
        best_step++;

      // compute new value
      if (formic::mpi::rank() == 0) {
        formic::xcopy(n, x0, 1, x1, 1);
        formic::xaxpy(n, step_sizes.at(step_i), p, 1, x1, 1);
      }
      worker.get_value(n, x1, step_energies.at(step_i));
      if (formic::mpi::rank() == 0) {
        formic::of << boost::format("for a step size of %16.8e") % step_sizes.at(step_i);
        formic::of << boost::format("      val = %20.12f") % step_energies.at(step_i) << std::endl << std::endl;
      }

      // if this was the best step size so far, save it
      if ( step_energies.at(step_i) < best_step_energy ) {
        best_step = step_i;
        best_step_size = step_sizes.at(step_i);
        best_step_energy = step_energies.at(step_i);
      }

      // if we have done two minima refinements, stop searching
      if (n_refine > 1)
        break;

      //// if the step sizes have become too close together, stop searching
      //const double step_thresh = 1.0e-6;
      //if ( best_step < step_sizes.size()-1 ) {
      //  if ( std::abs( ( step_sizes.at(best_step) - step_sizes.at(best_step+1) ) / step_sizes.at(best_step+1) ) < step_thresh )
      //    break;
      //}
      //if ( best_step > 0 ) {
      //  if ( std::abs( ( step_sizes.at(best_step) - step_sizes.at(best_step-1) ) / step_sizes.at(best_step) ) < step_thresh )
      //    break;
      //}

    }

    if (formic::mpi::rank() == 0) {

      // if a zero step was best, do not apply any update
      if ( best_step == 0 ) {

        formic::xcopy(n, x0, 1, x1, 1);

        formic::of << "The line search failed to find a step that lowers the energy.  No update was applied." << std::endl << std::endl;

      // otherwise, apply the update corresponding to the best step size
      } else {

        formic::xcopy(n, x0, 1, x1, 1);
        formic::xaxpy(n, best_step_size, p, 1, x1, 1);

        if (formic::mpi::rank() == 0)
          formic::of << boost::format("The line search's best step size was %12.6f with a value of %20.12f") % best_step_size % best_step_energy
                     << std::endl
                     << boost::format("The corresponding update has been applied.")
                     << std::endl
                     << std::endl;

      }

    }

    //// stop the timer
    //formic::stop_timer("lbgfs loose line search");

    // return the step size used
    step_arg = best_step_size;

  }

  template <class T> inline void lbfgs_shift_down(std::vector<T> & vec) {
      T temp = vec[0];
      for (int i = 0; i < vec.size()-1; i++)
        vec[i] = vec[i+1];
      vec[vec.size()-1] = temp;
  }

  //-------------------------------------------------------------------------------
  // formic::lbfgs -- This function performs an L-BFGS minimization of the
  //                  target function.
  //
  //  Inputs:
  //
  //    worker   -- object that computes the target function's value and gradient,
  //                and also deals with some other bookkeeping
  //    n        -- the dimension of the vectors x0 and x1
  //    x0       -- the initial guess.  upon exit, x0 is destroyed
  //
  //  Outputs:
  //
  //    x1       -- upon exit, the solution
  //
  //-------------------------------------------------------------------------------

  template <class W>
  inline void lbfgs(W & worker,
                    const int max_iter,
                    const int max_hist,
                    const size_t n,
                    double grad_thresh,
                    double value_thresh,
                    double * const x0,
                    double * const x1,
                    //double & step_size,
                    const bool print = true,
                    //const double c1 = 0.0001,
                    //const double c2 = 0.9,
                    const bool grad_by_ratio = false,
                    const int max_wolfe_iter = 10) {

    // get MPI info
    const int nproc = formic::mpi::size();
    const int myrank = formic::mpi::rank();

    // choose how many previous gradients to remember
    const int m = max_hist;

    // allocate storage for the gradient and position differences
    std::vector<double> grad_diff_vec((myrank == 0 ? size_t(n) * size_t(m) : size_t(1)), 0.0);
    std::vector<double>  pos_diff_vec((myrank == 0 ? size_t(n) * size_t(m) : size_t(1)), 0.0);
    std::vector<double *> grad_diffs(m);
    std::vector<double *> pos_diffs(m);
    if (myrank == 0) {
      for (size_t i = 0; i < size_t(m); i++) {
        grad_diffs[i] = &grad_diff_vec[i*n];
        pos_diffs[i] = &pos_diff_vec[i*n];
      }
    }

    // allocate gradients
    std::vector<double> grad_vec_0((myrank == 0 ? n : 1), 0.0);
    std::vector<double> grad_vec_1((myrank == 0 ? n : 1), 0.0);
    double * g0 = &grad_vec_0[0];
    double * g1 = &grad_vec_1[0];

    // allocate search direction
    std::vector<double> search_direction_vec((myrank == 0 ? n : 1), 0.0);
    double * p = &search_direction_vec[0];

    // allocate q vector
    std::vector<double> q_vec((myrank == 0 ? n : 1), 0.0);
    double * q = &q_vec[0];

    // allocate z vector
    std::vector<double> z_vec((myrank == 0 ? n : 1), 0.0);
    double * z = &z_vec[0];

    // allocate r and a vectors
    std::vector<double> r_vec((myrank == 0 ? m : 1 ), 0.0);
    std::vector<double> a_vec((myrank == 0 ? m : 1 ), 0.0);

    // prepare variables for the target function values
    double f0, f1;

    // compute the initial target function value
    worker.get_value(n, x0, f0);

    // compute the initial gradient
    worker.get_gradient(n, x0, g0);

    // print position and gradient
    worker.print(n, "initial position", x0);
    worker.print(n, "initial gradient", g0);

    // compute and print the norm of the initial gradient
    if (myrank == 0) {
      const double grad_norm = std::abs( std::sqrt( formic::xdotc(n, g0, 1, g0, 1) ) );
      if (print) {
        formic::of << boost::format(" function val = %20.12f") % f0 << std::endl;
        formic::of << boost::format("gradient norm = %20.12f") % grad_norm << std::endl << std::endl;
      }
      if (grad_by_ratio)
        grad_thresh *=  grad_norm;
    }

    // set the initial step size to try
    double step_size = 0.5;

    // iterate
    bool converged = false;
    int m_used = 0;
    for (int iter = 0; iter < max_iter; iter++) {

      //// set first step size to try
      //double step_size = 1.0;

      //// start a timer for the total iteration time
      //pcps::Stopwatch sw; sw.start();

      //// prepare a timer for sub-iteration timings
      //pcps::Stopwatch sw_sub;

      // print a greeting
      if (myrank == 0 && print)
        formic::of << boost::format("L-BFGS:  starting iteration %4i") % iter << std::endl << std::endl;

      //sw_sub.reset(); sw_sub.start();

      // compute the search direction on the root process ( p = - h * g0 )
      if (myrank == 0) {
        if (m > 0) {
          formic::xscal(n, 0.0, q, 1);
          formic::xaxpy(n, 1.0, g0, 1, q, 1);
          for (int i = m-1; i >= m - m_used; i--) {
            a_vec[i] = r_vec[i] * formic::xdotc(n, pos_diffs[i], 1, q, 1);
            formic::xaxpy(n, -a_vec[i], grad_diffs[i], 1, q, 1);
          }
          formic::xscal(n, 0.0, z, 1);
          formic::xaxpy(n, 1.0, q, 1, z, 1);  // assumes h_0 = identity
          for (int i = m - m_used; i < m; i++) {
            const double b = r_vec[i] * formic::xdotc(n, grad_diffs[i], 1, z, 1);
            formic::xaxpy(n, a_vec[i]-b, pos_diffs[i], 1, z, 1);
          }
          formic::xscal(n, 0.0, p, 1);
          formic::xaxpy(n, -1.0, z, 1, p, 1);
        } else {
          formic::xcopy(n, g0, 1, p, 1);
          formic::xscal(n, -1.0, p, 1);
        }
      }

      // print the search direction
      worker.print(n, "search direction", p);

      //sw_sub.stop();
      //if (myrank == 0 && print)
      //  formic::of << boost::format("prepare search direction took %.4f seconds") % sw_sub.elapsed_seconds() << std::endl;

      //// broadcast the search direction to all processes
      //pcps::bcast(comm, p, n, MPI::DOUBLE, 0);
      ////comm.Bcast((void *)p, n, MPI::DOUBLE, 0);

      //// perform a line search to find a new point x1 that satisfies the Wolfe conditions and the gradient g1 at that point
      //if (m == 0)
      //  pcps::sd_step_search(userinp, worker, n, c1, c2, step_size, p, x0, g0, f0, x1, g1, f1);
      //else 
      //  pcps::wolfe_search(userinp, worker, n, c1, c2, step_size, p, x0, g0, f0, x1, g1, f1, print, max_wolfe_iter);

      // perform a line search to find the updated position x1
      formic::quadratic_line_search(worker, n, f0, x0, p, x1, step_size);
      //formic::lbfgs_loose_line_search(worker, n, x0, p, x1, step_size);

      // if the step size was zero, we cannot descend
      formic::mpi::bcast(step_size);
      if ( myrank == 0 && print && step_size == 0.0 )
        formic::of << boost::format("computed search direction failed to lower the target function value, stopping lbfgs...") << std::endl << std::endl;
      if ( step_size == 0.0 )
        break;

      // compute new value
      worker.get_value(n, x1, f1);

      // compute new gradient
      worker.get_gradient(n, x1, g1);

      //// if necessary, increase the precision with which the values are calculated
      //worker.update_value_precision(f0, f1);

      // print the new position and gradient
      worker.print(n, "new position", x1);
      worker.print(n, "new gradient", g1);

      //sw_sub.reset(); sw_sub.start();

      // compute and print the norm of the new gradient
      const double grad_norm = ( myrank == 0 ? std::abs( std::sqrt( formic::xdotc(n, g1, 1, g1, 1) ) ) : 0.0 );
      if (myrank == 0 && print) {
        formic::of << boost::format(" function val = %20.12f") % f1 << std::endl;
        formic::of << boost::format("gradient norm = %20.12f") % grad_norm << std::endl << std::endl;
      }

      //if (myrank == 0 && grad_norm > userinp.max_grad_norm())
      //  pcps::dscal(n, 0.0, g1, 1);

      //sw_sub.stop();
      //if (myrank == 0 && print)
      //  formic::of << boost::format("gradient norm took %.4f seconds") % sw_sub.elapsed_seconds() << std::endl;

      //sw_sub.reset(); sw_sub.start();

      // if the gradient or energy difference is small enough, stop iterating
      converged = (grad_norm < grad_thresh) || (std::fabs(f1-f0) < value_thresh);
      //comm.Bcast((void *)&converged, sizeof(bool), MPI::CHAR, 0);
      formic::mpi::bcast(converged);
      if (converged) break;

      //sw_sub.stop();
      //if (myrank == 0 && print)
      //  formic::of << boost::format("broadcasting norm took %.4f seconds") % sw_sub.elapsed_seconds() << std::endl;

      // update the history of gradient differences, position differences, and r values
      if (myrank == 0 && m > 0) {

        lbfgs_shift_down(pos_diffs);
        formic::xscal(n, 0.0, pos_diffs[m-1], 1);
        formic::xaxpy(n, 1.0, x1, 1, pos_diffs[m-1], 1);
        formic::xaxpy(n,-1.0, x0, 1, pos_diffs[m-1], 1);
        //for (size_t j = 0; j < n; j++)
        //  pos_diffs[m-1][j] = x1[j] - x0[j];

        lbfgs_shift_down(grad_diffs);
        formic::xscal(n, 0.0, grad_diffs[m-1], 1);
        formic::xaxpy(n, 1.0, g1, 1, grad_diffs[m-1], 1);
        formic::xaxpy(n,-1.0, g0, 1, grad_diffs[m-1], 1);
        //for (size_t j = 0; j < n; j++)
        //  grad_diffs[m-1][j] = g1[j] - g0[j];

        lbfgs_shift_down(r_vec);
        r_vec[m-1] = 1.0 / formic::xdotc(n, grad_diffs[m-1], 1, pos_diffs[m-1], 1);

        m_used = std::min(m, m_used + 1);

      }

      //sw_sub.reset(); sw_sub.start();

      // perform any tasks requested to be done when a new position is found
      worker.process_new_point(n, x1);

      // compute new value now that we have processed the new point
      worker.get_value(n, x1, f1);

      // compute new gradient now that we have processed the new point
      worker.get_gradient(n, x1, g1);

      // move the new gradient, position, and function value into the old gradient, position, and function value
      if (myrank == 0) {
        f0 = f1;
        formic::xcopy(n, x1, 1, x0, 1); // x0 = x1
        formic::xcopy(n, g1, 1, g0, 1); // g0 = g1
      }

      //sw_sub.stop();
      //if (myrank == 0 && print)
      //  formic::of << boost::format("copying (x1, g1) to (x0, g0) took %.4f seconds") % sw_sub.elapsed_seconds() << std::endl;

      //// report the elapsed time for this iteration
      //sw.stop();
      //if (myrank == 0 && print)
      //  formic::of << boost::format("L-BFGS iteration took %.2f seconds") % sw.elapsed_seconds() << std::endl << std::endl;

    }

    // check that the method converged
    if (myrank == 0 && print) {
      if (converged)
        formic::of << boost::format("L-BFGS method converged!") << std::endl << std::endl;
      else
        formic::of << boost::format("L-BFGS method did not converge...") << std::endl << std::endl;
    }

  }

}

#endif
